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I. INTRODUCTION 



Utilization of quantum mechanical tunneling through 
classically localized states in annealing of glasses has 
opened up a new paradigm for solving hard optimiza- 
tion problems through adiabatic reduction of quantum 
fluctuations. This will be introduced and reviewed here. 

Consider the example of a ferromagnet consisting of 
N tiny interacting magnetic elements; the spins. For a 
macroscopic sample N is very large; of the order of Avo- 
gadro number. Assume that each spin can be in any of 
the two simple states: up or down. Also, the pairwise 
interactions between the spins are such that the energy 
of interaction (potential energy or PE) between any pair 
of spins is negative (smaller) if both the spins in the pair 
are in the same state and is positive (higher) if their 
states differ. As such, the collective energy of the TV- 
spin system (given by the Hamiltonian Ti.) is minimum 
when all the spins are aligned in the same direction; all 
up or all down, giving the full order. We call these two 
minimum energy configurations the ground states. The 
rest of the 2^ configurations are called excited states. 
The plot of the interaction energy of the whole system 
with respect to the configurations is called the potential 
energy-configuration landscape, or simply, the potential 
energy landscape (PEL). For a ferromagnet, this land- 
scape has a smooth double-valley structure (two mirror- 
symmetric valleys with the two degenerate ground states, 
all up and all down, at their respective bottoms). At zero- 
temperature the equilibrium state is the state of mini- 
mum potential energy, and the system stably resides at 
the bottom of any of the two valleys. At finite temper- 
ature, the thermal fiuctuations allow the system to visit 
higher energy configurations with some finite probabil- 
ity (given by the Boltzmann factor) and thus the system 
spends time in other part of the PEL also. The proba- 
bility that a system is found in a particular macroscopic 
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state depends not only on the energy of the state (as at 
zero temperature), but also on its entropy. The thermo- 
dynamic equilibrium state corresponds to the minimum 
of a thermodynamic potential called free energy F given 
by the difference between the energy of the state and 
the product of its entropy and the temperature. At zero- 
temperature, the minima of the free energy coincides with 
the minimum for energy and one gets the highest order 
(magnetization). As the temperature is increased, the 
contribution of entropy gets magnified and the minimum 
of the free energy is shifted more and more towards the 
states with lower and lower order or magnetization, un- 
til at (and beyond) some transition temperature Tc, the 
order disappears completely. For antiferromagnetic sys- 
tems, the spin-spin pair interactions are such that it is 
lower if the spins in the pair are in opposite states, and 
higher if their states are same. For antiferromagnets one 
can still define a sub-lattice order or magnetization and 
the PEL still has the double-well structure like in a fer- 
romagnet for short range interactions. The free energy 
and the order-disorder transition also shows identical be- 
havior as observed for a ferromagnet. 

In spin glasses, where different spin-spin interactions 
are randomly ferromagnetic or antiferromagnetic and 
frozen in time (quenched disorder), the PEL becomes ex- 
tremely rugged; various local and global minima trapped 
between potential energy barriers appears in the PEL. 
The ruggedness and the degeneracies in the minima 
comes from the effect of frustration or competing inter- 
actions between the spins; none of the spin states on a 
cluster or a placket is able to satisfy all the interactions 
in the cluster. The locally optimal state for the spins in 
the cluster is therefore degenerate and frustrated. 

Similar situations occur for multivariate optimization 
problems like the traveling salesman problem (TSP). 
Here a salesman has to visit A'' cities placed randomly 
on a plane (country). Of the N\/N distinct tours pass- 
ing through each city once, only few corresponds to the 
minimum (ground state) travel distance or travel cost. 
The rests corresponds to the higher costs (excited states). 
The cost function, when plotted against different tour 
configurations, gives a similar rugged landscape, equiva- 
lent to the PEL of a spin glass (henceforth we will use 
PEL to mean cost-configuration landscapes also). 

Obviously, an exhaustive search for the global mini- 
mum of a rugged PEL requires an exponentially large 
(or higher) number (in N) of searches (2^ or Nl or- 
der of searches for an A^-spin spin glass and an A-city 
TSP respectively). The computational effort or time for 
such searches are therefore generally not bounded by any 
polynomial in the problem size A. Alternatively, a grad- 
ual energy or cost dissipative dynamics (annealing), with 
Boltzmann like thermal fluctuations or some noise fac- 
tor (in order to get out of the the local minima) in the 
PEL may help reaching deep enough minima more eas- 
ily. This simulated thermal annealing scheme is consid- 
ered to be a very successful technique now. However, 
such a technique often fails if the barrier heights between 



the minima diverges (or become very high) , as in case of 
a spin glass, due to frequent trapping of the system in 
such local minima (glassy behavior). In case such bar- 
riers are very narrow, quantum-mechanical fluctuations 
(fluctuations in a quantum observable due to its non- 
commutativity with the Hamiltonian of the system) can 
help tunneling through them, thereby leading to success- 
ful quantum annealing (see. Fig. 1). 

We introduce here these ideas in details through suc- 
cessive steps. 

a) The physics of classical spin glasses has already con- 
tributed enormously to our knowledge of the landscape 
structure of the energy or the thermodynamic poten- 
tials and that of the unusually slow (glassy) dynamics 
of many-body systems in the presence of frustration and 
disorder. Mapping of computationally hard problems, 
like the traveling salesman problem etc, to classical spin 
glass models also helped understanding their complexity. 

b) The ground (and some low-lying) state structures of 
frustrated random systems in the presence of quantum 
fluctuations have also been studied in the context of 
quantum spin glasses. It has been shown that because 
of the possibility of tunneling through the barriers in 
the potential energy landscape, quantum fluctuations can 
help the dynamics to be 'more ergodic' than the dynam- 
ics induced by the classical fluctuations and thus help 
exploring the landscape much better. Ergodicity here 
means the loss of the memory of the initial state in course 
of evolution (weak ergodicity) and the convergence to a 
stationary distribution irrespective of the the initial state 
(strong ergodicity). The nature of these quantum phase 
transitions in such systems have also been extensively 
studied. These studies (Sec. IID) endow one with the 
knowledge of the phase diagram and the location of the 
quantum critical point or the phase boundary, which is 
crucial for choosing the proper quantum kinetic terms 
and the annealing path (Sec. IIIA and Sec. IIIB). 

c) The most natural connection between the paradigm 
of classical spin glasses and hard optimization problems 
comes through a widely used and well established op- 
timization technique, namely, simulated annealing algo- 
rithm as discussed earlier. The possibility of quantum 
tunneling through classically impenetrable barriers, as 
indicated from the studies of quantum spin glasses, nat- 
urally suggests an elegant and often more effective alter- 
native to simulated annealing. 

In quantum annealing, one has a classical Hamiltonian 
(or a multivariate cost function viewed as the same) to be 
optimized, to which, one adds a (non-commuting) quan- 
tum kinetic term and reduce it from a very high initial 
value to zero eventually. This reduction, when done com- 
pletely adiabatically, assures reaching of the ground state 
of the classical glass at the end, assuming that there is 
no crossing of energy levels with the ground state in the 
course of evolution, and provided that the starting state 
was the ground state of the initial Hamiltonian. To start 
with, the tunneling field is much higher than the inter- 
action term, so the ground state (a uniform superposi- 
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FIG. 1 While optimizing the cost function of a computationally 
hard problem (like the ground state energy of a spin glass or the 
minimum travel distance for a traveling salesman problem) , one has 
to get out of a shallower local minimum like the configuration C 
(spin configuration or travel route) , to reach a deeper minimum C' . 
This requires jumps or tunneling like fluctuations in the dynamics. 
Classically one has to jump over the energy or the cost barriers sep- 
arating them, while quantum mechanically one can tunnel through 
the same. If the barrier is high enough, thermal jump becomes 
very difficult. However, if the barrier is narrow enough, quantum 
tunneling often becomes quite easy. 

tion of all classical configurations) is trivially realizable. 
Simulations clearly demonstrate that quantum annealing 
can occasionally help reaching the ground state of a com- 
plex glassy system much faster than could be done using 
thermal annealing (discussed later in Sec. III). An exper- 
iment comparing the classical and quantum annealing for 
a spin glass also shows that the relaxations in course of 
quantum annealing are often much faster than those dur- 
ing the corresponding classical annealing, as discussed in 
Sec. HID. What makes quantum annealing fundamen- 
tally different from the classical annealing, is the non- 
local nature (Sec. Ill) and its higher tunneling ability 
(Sees. IID & IIIC). 

Quantum annealing thus permits a realization of ana- 
log quantum computation, which is an independent and 
powerful complement to digital quantum computation, 
where discrete unitary transformations are implemented 
through quantum logic gates. 

II. OPTIMIZATION AND ANNEALING 

A. Combinatorial Optimization Problems 

The occurrence of multivariate optimization problems is 
ubiquitous in our life, wherever one has to choose the 
best bargain from a host of available options that de- 
pend on many independent factors. In many cases, such 
a task can be cast as a problem of minimizing a given 
cost or energy function T-L{Si,S2, ■■■Sn) with respect to 
N variables Si, S2, ■■■Sn (sometimes subject to some con- 
straints). The task is to find a set of values for these vari- 
ables (a configuration) for which the function T-C{{Si}) 
has the minimum value (cf. Fig. 1). In many important 
optimization problems, the set of feasible configurations 
from which an optimum is to be chosen is a finite set (for 
finite N). In such a case, we say that the problem is com- 



binatorial in nature. If the variables Si are discrete and 
each takes up a finite number of values, then the problem 
is clearly a combinatorial one. Moreover, certain prob- 
lems with continuous variables (like linear programming 
problem) can also be reduced to combinatorial problems 
(Papadimitriou et al 1998). Here we focus on this type 
of optimization problem, and assume that we have to 
minimize H{{Si}) with respect to the discrete set of the 
variables Si. 

An optimization problem is said to belong to the class 
P (P for Polynomial), if it can be solved in polynomial 
time (i.e., the evaluation time goes like some polynomial 
in N) using polynomially (in N, again) bound resources 
(computer space, processors etc). Existence of such a 
polynomial bound on the evaluation time is somehow in- 
terpreted as the "easiness" of the problem. However, 
many important optimization problems seem to fall out- 
side this class, like, the famous traveling salesman prob- 
lem (see Sec. HC.2). 

There is another important class of problems which can 
be solved in polynomial time by non-deterministic ma- 
chines. This class is the famous NP (Non-deterministic 
Polynomial) class (Garey and Johnson 1979). P is in- 
cluded completely in the NP class, since a deterministic 
Turing machine is just a special case of non-deterministic 
Turing machines. Unlike a deterministic machine, which 
takes a specific step deterministically at each instant 
(and hence follows a single computational path), a non- 
deterministic machine has a host of different 'allowed' 
steps at its disposal at every instant. At each instant 
it explores all the 'allowed' steps and if any of them 
leads to the goal, the job is considered to be done. Thus 
it explores in parallel many paths (whose number goes 
roughly exponentially with time) and checks if any one 
of them reaches the goal. 

Among the NP problems, there are certain problems 
(known as NP-complete problems) which are such that 
any NP problem can be "reduced" to them using a poly- 
nomial algorithm. The famous 3-SAT problem (see Sec. 
IIIA.3) is a representative of the class. This roughly 
means that if one has a routine to solve an NP-complete 
problem of size N then using that routine one can solve 
any NP problem at the cost of an extra overhead in time 
that goes only polynomially with A^. The problems in 
this class are considered to be hard, since so far no one 
can simulate a general nondeterministic machine by a de- 
terministic Turing machine (or any sequential computer 
with polynomially bound resources) without an exponen- 
tial growth of execution time. In fact it is largely believed 
(though not proved yet) that it is impossible to do so (i.e., 
Py^NP) in principle. However, assuming this to be true, 
one can show that there are indeed problems in NP class 
that are neither NP-complete nor P (Garey and Johnson 
1979). 
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B. Statistical IVIechanics of the Optimization Problems and 
Thermal Annealing 

There arc some excellent deterministic algorithms for 
solving certain optimization problems exactly (Papadim- 
itriou and Steiglitz 1998, Hartmann and Rieger, 2002). 
These algorithms are, however, quite small in number 
and are strictly problem specific. For NP or harder prob- 
lems, only approximate results can be found using these 
algorithms in polynomial time. These approximate algo- 
rithms too are also strictly problem specific in the sense 
that if one can solve a certain NP-completc problem up 
to a certain approximation using some polynomial algo- 
rithm, then that does not ensure that one can solve all 
other NP problems using it up to the said approximation 
in polynomial time. 

Exact algorithms being scarce, one has to go for heuris- 
tics algorithms, which are algorithms based on certain 
intuitive moves, without any guarantee on either the 
accuracy or the run time for the worst case instance. 
However, these algorithms arc generally easy to formu- 
late and are quite effective in solving most instances of 
a the intended problems. A general approach towards 
formulating such approximate heuristics may be based 
on stochastic (randomized) iterative improvements. The 
most preliminary one in this family is the local mini- 
mization algorithm. In this algorithm one starts with a 
random configuration Co and makes some local changes 
in the configuration following some prescription (stochas- 
tic or deterministic) to generate a new configuration Ci 
and calculates the corresponding change in the cost. If 
the cost is lowered by the change, then the new config- 
uration Ci is adopted. Otherwise the old configuration 
is retained. Then in the next step a new local change 
is attempted again, and so on. This reduces the cost 
steadily until a configuration is reached which minimizes 
the cost locally. This means that no further lowering of 
cost is possible by changing this configuration using any 
of the prescribed local moves. The algorithm essentially 
stops there. But generally in most optimization problems 
(such as in spin glasses), there occur many local minima 
in the cost-configuration landscape and they are mostly 
far above the global minimum (see Fig. 1). It is likely 
that the algorithm therefore gets stuck in one of them and 
ends up with a poor approximation. One can then start 
afresh with some new initial configuration and end up 
with another local minimum. After repeating this for sev- 
eral times, each time with a new initial configuration, one 
may choose the best result from them. But much better 
idea would be to get somehow out of shallow local min- 
ima. One can introduce some fluctuations or noise in the 
process so that the movement is not always towards lower 
energy configurations, but there is also a finite probabil- 
ity to go to higher energy configurations (the higher the 
final energy, the lower the probability to move to that), 
and consequently there appear chances to get out of the 
shallow local minima. Initially, strong fiuctuations are 
adopted (i.e., the probability to go to higher energy con- 



figurations is relatively high) and slowly fiuctuations are 
reduced until finally they are tuned off completely. In 
the mean time the system gets a fair opportunity to ex- 
plore the landscape more exhaustively and settle into a 
reasonably deep cost or energy minimum. Kirkpartick 
et al (1983) suggested an elegant way: A fluctuation is 
implemented by introducing an "artificial" temperature 
T into the problem such that the transition probability 
from a configuration Ci to a configuration C/ is given 
by min{l,exp — [Aj//T]}, where Aj/ = Ef — Ei, with 
Ek denoting the cost or energy of the configuration Cfc. 
A corresponding Monte Carlo dynamics is defined, say, 
based on detailed balance, and the thermal relaxation of 
the system is simulated. In course of simulation, the noise 
factor T is reduced slowly from a very high initial value 
to zero following some annealing schedule. At the end of 
the simulation one is expected to end up with a config- 
uration whose cost is a reasonable approximation of the 
globally minimum one. If the temperature is decreased 
slow enough, say. 



T{t)>N/logt, (1) 



where t denotes the cooling time and A'' the system size, 
then the global minimum is attained with certainty in the 
limit t oo (Geman and Geman 1984). Even within a 
finite time and with a faster cooling rate, one can achieve 
a reasonably good approximation (a crystal with only few 
defects) in practice. This simulated annealing method is 
now being used extensively by engineers for devising real- 
life optimization algorithms. We will refer to this as the 
classical annealing (CA), to distinguish it from quantum 
annealing (QA) which employs quantum fluctuations. It 
is important to note that though in this type of stochastic 
algorithms the system has many different steps with their 
respective probabilities at its disposal, it finally takes up 
a single one, say by tossing coins, and thus finally follow 
a single (stochastically selected) path. Hence it is not 
equivalent to a non-deterministic machine, where all the 
allowed paths are checked in parallel at every time-step. 

As has been mentioned already, many combinatorial 
optimization problems can be cast into the problem of 
finding the ground state of some classical (spin glass like) 
Hamiltonian H{{Si}). One can therefore analyze the 
problem by using statistical mechanics so as to to ap- 
ply physical techniques like simulated annealing. If one 
naively takes the number of variables A'' as the size, then 
the entropy and the energy are often found to scale dif- 
ferently with N and the application of standard thermo- 
dynamic arguments become difficult. One needs to scale 
temperature and some other quantities properly with N 
so that one can talk in terms of the concepts like free en- 
ergy minimization etc. Moreover, the constraints present 
in the problems are often very difficult to take into ac- 
count. 
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C. Spin Glass and Optimization 

1. Finding The Ground States of the Classical Spin Glasses 

As has been mentioned already, the difficulty faced by 
a physically motivated optimization heuristic (one that 
follows physical relaxation dynamics, classical or quan- 
tum, to search for the solution), in finding the solution 
of a hard optimization problem, is similar to that faced 
by a glassy system in reaching its ground state. In fact, 
finding the ground states of the spin glasses is an impor- 
tant class of combinatorial optimization problem, which 
includes an NP-completc problem (Barahona 1982), and 
many other apparently different ones (like the traveling 
salesman problem) can be recast in this form. Hence here 
we discuss very briefly the nature of the spin glass phase 
and the difficulty in reaching its ground state. 

The interaction energy of a typically random and frus- 
trated Ising spin glass (Binder and Young 1986, Nishi- 
mori 2001, Dotsenko 2001) may be represented by a 
Hamiltonian of the form 

N 

H = - ^JijSiSj, (2) 

i>j 

where Si denote the Ising spins and the interactions 
between them. These JjjS here are quenched variables 
which vary randomly both in sign and magnitude follow- 
ing some distribution p(Jij). The typical distributions 
are (zero mean) Gaussian distribution of positive and 
negative values 

p(Ji,)=Aexp|^^j, (3) 

and binary distribution 

piJij) = p5{.hj - J) + (1 - p)5{J,j + J), (4) 

with probability p of having a + J bond, and (1 — p) of 
having a — J bond. Two well studied models are the 
SK model (introduced by Sherrington and Kirkpatrick 
1975) and the EA model (introduced by Edward and 
Anderson 1975). In the SK model the interactions are 
infinite-ranged and for the sake of extensively (for rules 
of equilibrium thermodynamics to be applicable, energy 
should be proportional to the volume or its equivalent 
that defines the system-size) one has to scale J 1/ y/N, 
while in the EA model, the interactions are between the 
nearest neighbors only. For both of them, however, p{Jij) 
is Gaussian; A = (A^/27rJ^)^/^ for normalization in the 
SK model. 

The freezing (below Tc) is characterized by some non- 
zero value of the thermal average of the magnetization 
at each site (local ordering). However, since the interac- 
tions are random and competing, the spatial average of 
single-site magnetization (below Tc) is zero. Above Tc, 



both the spatial and the temporal averages of the single- 
site magnetization of course vanishes. A relevant order 
parameter for this freezing is therefore 
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9=7V^<^^^T, (5) 

(overhead bar denoting the average over disorders (the 
distribution of Jy ) and (...)t denotes the thermal aver- 
age). Here q ior T < Tc, while g = for T > T^. 
As will be seen in the following, the existence of a unique 
order parameter q will indicate ergodicity. 

In the spin glass phase (T < Tc), the whole free energy 
landscape gets divided (cf. Fig. 1) into many valleys 
(local minima of free energy) separated by very high free 
energy barriers. Thus the system, once trapped in a val- 
ley, remains there for a very long time. The spins of such 
a confined system are allowed to explore only a very re- 
stricted (and correlated) part of the configuration space 
and thus makes them "freeze" to have a magnetization 
that characterizes the state (valley) locally. 

So far, two competing pictures continue to represent 
the physics of the spin glasses. The mean-field picture 
of replica symmetry-breaking is valid for infinite-ranged 
spin glass systems like the SK spin glass. In this pic- 
ture, below the glass transition temperature Tc, the bar- 
riers separating the valleys in the free energy landscape 
actually diverge (in the limit N — > oo), giving rise to 
a diverging timescale for the confinement of the system 
in any such valley once the system gets there somehow. 
This means there is a loss of ergodicity in the thermal 
dynamics of the system at T < Tg. Thus one needs a 
distribution P{q) of order parameters, instead of a sin- 
gle order parameter to characterize the whole landscape, 
as emerges naturally from the replica symmetry-breaking 
ansatz of Parisi (1980). To be a bit more quantitative, let 
us imagine that two identical replicas (having exactly the 
same set J^ s) of a spin glass sample are allowed to relax 
thermally below Tc, starting from two different random 
(paramagnetic) initial states. Then these two replicas 
(labeled by /x and i^, say) will settle in two different val- 
leys, each being characterized by a local value of the order 
parameter and the corresponding overlap parameters q^v 
having a sample-specific distribution 

Pjil) =^e—^S{q-qi^^); 

P{q) = ll[dJijp{Jij)Pj{q). (6) 

Here the subscript J denotes a particular sample with a 
given realization of quenched random interactions (Jijs) 
between the spins, and finally averaging Pj{q) over dis- 
order distribution p{J) in (2) or (3) one gets P{q). Phys- 
ically, P{q) gives the probability distribution for the two 
pure states to have an overlap q, assuming that the proba- 
bility of reaching any pure state p starting from a random 
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(high temperature) state is proportional to the thermo- 
dynamic weight exp {— -F),} of the state fi. The other pic- 
ture is due to the droplet model of short-range spin glass 
(Bray and Moore 1984, Fisher and Huse 1986), where 
there is no divergence in the typical free energy barrier 
height and the relevant timescale is taken to be that of 
crossing the free energy barrier of formation of a typi- 
cal droplet of same (all up or all down) spins. Based on 
certain scaling ansatz, this picture leads to a logarith- 
mically decaying (with time) self-correlation function for 
the spins below the freezing temperature T^. 

The issue of the validity of the mean-field picture (of 
the replica symmetry-breaking) in the context of the real- 
life spin glasses, where the interactions are essentially 
short-range, is far from settled (see e.g., Moore et al 1998, 
Marinari et al 1998, Krzakala et al 2001 and Gaviro et 
al 2006). However, the effective Hamiltonian (cost func- 
tion) for many other optimization problems may contain 
long-range interactions and may even show the replica 
symmetry-breaking behavior shown in the Graph Parti- 
tioning Problem (Fu and Anderson 1986). Of course, no 
result of QA for such system (for which replica symmetry- 
breaking is shown explicitly) has been reported yet. The 
successes of QA reported so far are mostly for short-range 
systems. Thus, the scope of quantum annealing in those 
long-range systems still stands out as an interesting open 
question. 



2. The Traveling Salesman Problem 

In the traveling salesman problem (TSP), there arc N 
cities placed randomly in a country having a definite met- 
ric to calculate the inter-city distances. A salesman has 
to make a tour to cover every city and finally come back 
to the starting point. The problem is to find a tour of 
minimum length. An instance of the problem is given by 
a set {dij-, i,j = 1, A^}, where dij indicates the distance 
between the i-th and the j-th city, or equivalently, the 
cost for going from the former to the later. We mainly 
focus on the results of symmetric case, where dij = dji. 
The problem can be cast into the form where one min- 
imizes an Ising Hamiltonian imdcr some constraints, as 
shown below. A tour can be represented by an A x A^ 
matrix T with elements either or 1. In a given tour, 
if the city j is visited immediately after visiting city i, 
then Tij = 1, or else = 0. Generally an additional 
constraint is imposed that one city has to be visited once 
and only once in a tour. Any valid tour with the above 
restriction may be represented by a T matrix whose each 
row and each column has one and only one element equal 
to 1 and rest all are Os. For symmetric a metric, a 
tour and its reverse tour have the same length, and it is 
more convenient to work with an undirected tour matrix 
U = ^(T -|- T), where T, the transpose of T, represents 
the reverse of the tour given by T. Clearly, U must be a 
symmetric matrix having two and only two distinct en- 
tries equal to 1 in every row and every column, no two 



rows being identical, and so is not any two columns. In 
terms of UijS, the length of a tour can be represented by 

1 ^ 

E'^^^U^.- (7) 

One can rewrite the above Hamiltonian in terms of Ising 
spins SijS as 

where Sij = 2Uij — 1 are the Ising spins. The Hamilto- 
nian is similar to that of a non-interacting Ising spins on 
an N X N lattice, with random fields dij on the lattice 
points {i,j}- The frustration is introduced by the global 
constraints on the spin configurations in order to conform 
with the structure of the matrix U discussed above. The 
problem is to find the ground state of the Hamiltonian 
subject to these constraints. There are Ising spins, 
which can assume 2^ configurations in absence of any 
constraint, but the constraint here reduces the number 
of valid configurations to that of the number of distinct 
tours, which is {N\)/2N. 

Mainly two distinct classes of TSP arc studied: one 
with an Euclidean dij in finite dimension (where dij are 
strongly correlated through triangle inequalities, which 
means, for any three cities A, B and C, the sum of any 
two of the side AB, BC and CA must be greater than 
the remaining one), and the other with random dij in 
infinite dimension. 

In the first case, N cities are uniformly distributed 
within a hypercube in a d-dimensional Euclidean space. 
Finding a good approximation for large N is easier for 
this case, since the problem is finite-ranged. Here a d- 
dimensional neighborhood is defined for each city, and 
the problem can be solved by dividing the whole hyper- 
cube into a number of smaller pieces and then searching 
for the least path within each smaller part and joining 
them back. The correction to the true least path will be 
due to the unoptimized connections across the bound- 
aries of the subdivisions. For a suitably made division 
(not too small), this correction will be of the order of 
the surface-to-volume ratio of each division, and thus 
will tend to zero in the N ^ 00 limit. This method, 
known as "divide and conquer" , forms a reasonable strat- 
egy for solving approximately such finite-range optimiza- 
tion problems (including finite-range spin glasses) in gen- 
eral. In the second case, d^ s are assigned completely ran- 
domly, with no geometric (e.g., Euclidean) correlation be- 
tween them. The problem in this case becomes more like 
a long-range spin glass. A self avoiding walk represen- 
tation of the problem was made using an m-component 
vector field, and the replica analysis was done (Mezard et 
al 1987) for finite temperature, assuming replica symme- 
try ansatz to hold. Moreover, true breaking of ergodicity 
may occur only in infinite systems, not in any finite in- 
stance of the problem. The results, when extrapolated to 
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zcro-tcmpcraturc, do not disagree much with the numer- 
ical results (Mezard et al 1987). The stability of a replica 
symmetric sohition has not yet been proved for low tem- 
perature region. However, numerical results of thermal 
annealing for instances of size A'' = 60 to A'' = 160 yielded 
many near optimal tours, and the corresponding over- 
lap analysis shows a sharply peaked distribution, whose 
width decreases steadily with increase in N. This indi- 
cates the existence of a replica symmetric phase for the 
system (Mezard et al 1987). 

An analytical bound on the average (Normalized by 
N^^"^) value of optimal path-length per city (Q) calcu- 
lated for TSP on 2-dimensional Euclidean plane has been 
found to be 5/8 < < 0.92 (Bearwood et al 1959). Care- 
ful scaling analysis of the numerical results obtained so 
far indicates the lower bound to be close to 0.72 (Percus 
and Martin 1996, Chakraborti and Chakrabarti 2000). 

Simulated (thermal) annealing of a Euclidean TSP on 
a square having length N^^"^ (so as to render the average 
nearest neighbor distance independent of N) has been re- 
ported (Kirkpatrick et al 1983). In this choice of length 
unit, the optimal tour length per step (n) becomes inde- 
pendent of N for large N. Thermal annealing rendered 
n < 0.95 for N up to 6000 cities. This is much better 
than what is obtained by the so called "greedy heuris- 
tics" , (where being at some city in a step, one moves to 
the nearest city not in the tour in the next step) for which 
O ^ 1.12 on average. Later we will see that (Sec. IIIA) 
that quantum annealing can do even better than thermal 
annealing in context of random TSP. 

To summarize, when cast in some energy minimiza- 
tion problem, the combinatorial optimization problems 
may exhibit glassy behavior during thermal annealing. 
Even replica symmetry-breaking behavior may be ob- 
served (like in case of Graph Partitioning Problem, see Fu 
and Anderson 1986), since the underlying Hamiltonian 
need not be short-range, and the constraints can bring 
frustrations into the problem. One can intuitively con- 
clude that thermal annealing or other heuristics would 
not be able to solve such problems easily to a good ap- 
proximation within reasonable time. Moreover, practi- 
cally nothing can be said about the time required to 
solve the worst case instance exactly. Specifically, in 
some cases, where good solutions are thermodynamically 
very insignificant in number and there is no monotonic 
gradient towards them, the entropy might make a clas- 
sical search exponentially difficult, though the landscape 
might still remain completely ergodic. Later we will sec 
quantum searches can bring about spectacular improve- 
ments in some such cases (see Sec. IIIB and Fig. 4). 



D. Quantum Spin Glasses and Annealing 

In QA one adds a kinetic (tunneling) term to the inter- 
action part of the classical glass Hamiltonian. The ob- 
ject that results in, is called a quantum spin glass. The 
knowledge of the phase-diagram of a quantum spin glass 



is crucially important for its annealing, as it gives an 
idea of the location of the quantum critical points on the 
phase diagram, and thus offers a guideline to choose the 
proper kinetic terms (that maintain a sizable gap) and 
the suitable annealing paths (see Sec. IIIA and IIIB). 

In quantum spin glasses (Chakrabarti 1981, Ishii and 
Yamamoto 1985, Bhatt 1998, Rieger 2005, Sachdev 1999, 
Ye 1993), the order-disorder transition (i.e., from the 
frozen phase to the high kinetic energy phase, termed 
the para phase) can be driven both by thermal fluctua- 
tions as well as by quantum fluctuations. Quantum spin 
glasses can be of two types: vector spin glasses, where 
the quantum fluctuations cannot be adjusted by chang- 
ing some laboratory field, and the other, a classical spin 
glass perturbed by some quantum tunneling term, where 
the quantum fluctuations are controlled through, say, a 
transverse laboratory field. 

The amount of the quantum fluctuations being ad- 
justable, this transverse Ising spin glass (TISG) model 
is perhaps the simplest model in which the quantum ef- 
fects in a random system can be and have been studied 
extensively and systematically (Chakrabarti et al 1996). 
Here we will focus only on TISG, since the reduction of 
the quantum fluctuations is the key feature required for 
quantum annealing. 

The interest in the zero-temperature quantum spin 
glass phases in the TISG models have been have been 
complimented all along by the experimental studies 
in several systems which have been established to be 
represented accurately by transverse field Ising model 
(TIM). Recent discovery of the compoimd material 
LiHoa;Yi_a;F4 with the magnetic Ho ion concentration 
X = 0.167 (Aepph and Rosenbaum 2005, Wu et al 1991, 
Wu et al 1993a; see also Brooke et al 2001, Silevitch et al 
2007) representing accurately a random long-range trans- 
verse Ising system, has led to renewed interest. Here, the 
strong spin-orbit coupling between the spins and the host 
crystals restricts the effective "Ising" spins to align either 
parallel or anti-parallel to the specific crystal axis. An ap- 
plied magnetic field, transverse to the preferred axis, flips 
the "Ising" spins. This feature, together with the ran- 
domness in the spin-spin interaction, makes it a unique 
TISG-like system. Most interestingly, it has been shown 
that in spite of the presence of all the three ingredients 
- frustrations, randomness and the long-range (dipolar) 
interactions, that are necessary for the formation of a 
spin glass, the spin glass phase of LiHo2;Yi_a;F4 is ac- 
tually destroyed by any finite transverse field (Schechter 
and Lafiorencie 2006). This indicates the effectiveness of 
quantum tunneling in the exploration of a rugged PEL 
with formidable potential energy barriers. The TISG 
model described here, is given by the Hamiltonian 

JV N 

n = -Y,JioS!S]-TY,sf, (9) 

i>i i 

where V denotes the tunneling strength at each site and 
JijS are distributed randomly following the distribution 
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p{Jij) given by (3) or (4). Generally, we denote the 
strength of the quantum kinetic term by T. 

The unique interest in such quantum spin glass sys- 
tem comes from the possibility of much faster crossing of 
the high barriers occurring in the potential energy land- 
scapes of the classical spin glasses by means of quantum 
tunneling induced by the transverse field, compared to 
that done thermally by scaling such barriers with the 
temperature. 

The phase transitions in quantum spin glasses can be 
driven both by thermal and quantum fluctuations as 
mentioned before and the equilibrium phase diagrams 
also indicate how the optimized solution (in SG phase) 
can be obtained either by tuning of the temperature T or 
the tunneling field F, or by both. We will show later (in 
context of quantum annealing) that reaching the phase 
by tuning T may often be more advantageous (faster) 
than that by tuning T. 

The short-range version of this TISG model was first 
studied by Chakrabarti (1981), and the long-range ver- 
sion, discussed here, was first studied by Ishii and Ya- 
mamoto (1985). Several analytical studies have been 
made to obtain the phase diagram of the transverse Ising 
SK model (Miller and Huse 1993). The problem of an 
SK glass in a transverse field becomes a nontrivial one 
due to the presence of non-commuting spin operators in 
the Hamiltonian. This leads to a dynamical frequency- 
dependent self- interaction for the spins. 

One can study an effective-spin Hamiltonian for the 
SK model in a transverse field within the mean field 
framework very easily. The spin glass order parameter 
in a classical SK model is given by a random 'mean field' 
h{r) having Gaussian distribution (see Binder and Young 
1986) 

/ + 00 
dre-'''/2 tanh^ (h''{r)/T); ft,^(r) = Jy^r+h' 
-OO 

(10) 

where denotes the external field (in z direction), the 
mean-field h{r) being also in the same direction. In the 
presence of the transverse field, as in (9), h{r) has com- 
ponents both in z and x directions 

h{r) = -h''{r)z - Fx; h{r) = ^/h^{ry +T^, (11) 

and one replaces the ordering term tanh^ {h{r)/T) in (10) 
by its component Wh"" {r)\/\h{r)\f tanh^ {\h{r)\/T) in the 
2:-direction. Putting =Q and Q one gets the phase 
boundary equation as (see Chakrabarti et al 1996) 

^=tanh(^^). (12) 

This gives {T = 0) = J = Tc {T = 0) and a 
phase diagram qualitatively similar to the experimental 
one shown in Fig. 2. 

Several Monte Garlo studies have been performed for 

the SK spin glass model in transverse field applying the 
Suzuki- Trotter formalism (see Appendix 1), mapping a 



rf-dimensional quantum Hamiltonian to an effective d+1 
dimensional anisotropic classical Hamiltonian (see also 
Hatano and Suzuki 2005). The partition function gives 
the effective classical Hamiltonian in the Mth Trotter 
approximation as 

N M N M 

= ^ ^ KijSikSjk — ^ ^ KSikSik+i, (13) 

i>j k i k 

with 

= ^=^l-oth(^), (14) 

where Si^ denotes the Ising spin defined on the lattice 
site (i, fc), i denoting the position in the in the original 
SK model and k denoting the position in the additional 
Trotter dimension. Although the equivalence between 
classical and the quantum model holds exactly in the 
limit M ^ OO, one can always make an optimum choice 
for M. The equivalent classical Hamiltonian has been 
studied using standard Monte Carlo technique. The nu- 
merical estimates of the phase diagram etc. are reviewed 
in details in (Bhatt 1998, Rieger 2005). Ray et al (1989) 
took F << J and their results indeed indicated a sharp 
lowering of Tc{T). Such sharp fall of Tc(F) with large 
r is obtained in almost all the theoretical studies of the 
phase diagram of the model (Miller and Huse 1993; Ye et 
al 1993; see also Bhatt 1998 and Rieger 2005). Quantum 
Monte Carlo (Alverez and Ritort 1996) as well as real- 
time Schrodinger evolution (the true dynamics given by 
the time-dependent Schrodinger equation) studies of SK 
spin glass in transverse field were made (Lancaster and 
Ritort, 1997). 

In the Hamiltonian for the EA spin glass in presence of 
a transverse field, given by (9), the random interactions 
are restricted among the nearest neighbors and satisfy 
a Gaussian distribution with zero mean and variance J, 
as given by Eq. (3). Here, the variation of correlations 
in the equivalent {d+1) dimensional classical model fit- 
ted very well (Guo et al 1994) with the scaling fit with 
a unique order parameter and a critical interval corre- 
sponding to a phase diagram whose features are similar 
to those discussed above (see also Chakrabarti et al 1996, 
Bhatt 1998 and Rieger 2005). 

As discussed earlier in this section, LiHoj;Yi_j.F4 with 
X = 0.167 provides a spin glass system, for which the 
external magnetic field transverse to the preferred axis 
scales like the square root of the tunneling field F in 
(9). With increasing transverse field, the glass transition 
temperature decrease monotonically, as shown in Fig. 2. 
A quantum tunneling term allows for overlap between 
two classically localized states and the dynamics near 
the ground states of such glasses show better ergodicity 
properties. In order to investigate this aspect of quan- 
tum spin glasses, one can study the overlap distribution 
function P{q) given by Eq. (6). 
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FIG. 2 Phase diagram of LiHoo.i67Yo.833F4 according to the 

dynamical measurements (filled circles) and of the nonlinear sus- 
ceptibility measurements (open circles). Filled squares indicate the 
freezing boundary obtained from AC susceptibility measurements 
(taken from Wu et al 1993b). 

If the crgodicity is recovered, at least for a part of the 
phase diagram, the above function should tend to a delta 
function form, peaking at some finite value of the order- 
parameter q in thermodynamic limit. In the para-phase, 
of course, the distribution becomes a delta function at 
g = for the infinite system. In spite of several investi- 
gations ( see e.g., Ray et al 1989, Thirumalai and Kirk- 
patrick 1989, Goldschmidt and Lai 1990, Chakrabarti et 
al 1996, Kim and Kim 2002) the point of replica symme- 
try restoration in spin glasses by quantum fluctuation is 
not settled yet. However, slow withdrawal (sec Eq. 16 
for the characteristic slowness) of the tunneling field in 
these quantum spin glasses can help annealing the sys- 
tem close to the ground state of the classical spin glass 
eventually, as has been described in the next section. 

III. QUANTUM ANNEALING 

In the previous sections we have seen how thermal fluctu- 
ations can be utilized to devise fast heuristics to find an 
approximate ground state of a glassy system, or equiva- 
lently, a near-optimal solution to combinatorial problem, 
whose cost-configuration landscape has glassy behavior 
due to the occurrence of many local minima. There 
are two aspects of an optimization problem which might 
render thermal annealing to be a very ineffective one. 
First, in a glassy landscape, there may exist very high 
cost/energy barriers around local minima which does not 
correspond to a reasonably low cost (sec Fig. 1). In the 
case of infinite-range problems, these barriers might be 
proportional to the system-size iV, and thus diverge in 
thermodynamic limit. Thus there might occur many un- 
satisfactory local minima, any of which can trap the sys- 
tem for very long time (which actually diverges in ther- 
modynamic limit for infinite-range systems) in course of 
annealing. The second problem is the entropy itself. The 
number of configurations grow very fast with the number 
of variables (roughly exponentially; n Ising spin can be in 
2" configurations) and a classical system can only assume 
one configuration at a time and unless there is a gradi- 
ent that broadly guides the system towards the global 



minimum form any point in the configuration space, the 
search has to involve visiting any substantial fraction of 
the configurations. Thus a PEL without a guiding gra- 
dient poses a problem which is clearly of an exponen- 
tial or higher order in complexity (depending how the 
size of the configuration space scales with the system- 
size N), and the CA algorithms can do no better than 
a random search algorithm. This is the case for golf- 
course type potential-energy landscape (where there is a 
sharp potential minimum on completely fiat PEL) (see 
Sec. IIIB, Fig. 4). One can imagine that quantum me- 
chanics might have some solutions to both these prob- 
lems, at least, to some extent. This is because quan- 
tum mechanics can introduce classically unlikely tunnel- 
ing paths even through very high barriers if they are nar- 
row enough (Ray et al 1989; see also Apolloni et al 1989, 
1990). This can solve the ergodicity problem to some ex- 
tent, as discussed earlier. Even in places where ergodicity 
breaking does not take place in true sense, once the en- 
ergy landscape contains high enough barriers (specially 
for infinite-ranged quenched interactions), quantum tun- 
neling may provide much faster relaxation to the ground 
state (Santoro et al 2002, Martofiak et al 2002; see also 
Santoro and Tosatti 2007, Somma et al 2007). In addi- 
tion, a quantum mechanical wave function can delocalize 
over the whole configuration space (i.e., over the poten- 
tial energy landscape) if the kinetic energy term is high 
enough. Thus it can actually "see" the whole landscape 
simultaneously at some stage of annealing. These two 
aspects can be expected to improve the search process 
when employed properly. In fact such improvements can 
indeed be achieved in certain situations, though quan- 
tum mechanics is not a panacea for all such diseases as 
ergodicity breaking, spin glass behavior etc, and certainly 
has its own inherent limitations. What intrigues one is 
the fact that the limitations due to quantum nature of 
an algorithm are inherently different from those faced 
by its classical counterpart, and thus it is not yet clear 
in general which wins when. Here we discuss some re- 
sults regarding the quantum heuristics and some of their 
general aspects understood so far. For more detailed re- 
views of the subject wc refer to the articles in Das and 
Chakrabarti (2005) and Santoro and Tosatti (2006). 

Some basic aspects of QA can be understood from the 
simple case of QA in the context of a double-well po- 
tential (Stella 2005, Battagha et al 2005). Typically a 
particle in a double-well consisting of a shallower but 
wider well and a deeper but narrower well is annealed 
(it is likely that the deeper well, i.e., the target state, 
is narrower, otherwise the searching becomes easier even 
classically). The kinc^tic cinergy (inverse mass) is tuned 
from a very high value to zero linearly within a time r. 
For a very high value of initial kinetic energy, the wave 
function, which is the grotmd state, is delocalized more 
or less over the whole double-well. As kinetic energy is 
reduced but still quite high, the ground state corresponds 
to a more pronounced peak on the shallower minimum, 
since it is wider. This is because at this stage, to obtain 
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the minimal (ground state) energy, it is more effective to 
minimize the kinetic energy by localization over a wider 
space, rather than minimizing the potential energy by 
localizing in the deeper well. However, as kinetic en- 
ergy is tuned down further, the potential energy term 
becomes dominating, and the ground state has a taller 
peak around the deeper minimum. The evolving wave 
function can roughly follow this ground state structure 
all the way and finally settle to the deeper minimum if 
the annealing time t is greater than some r^. When 
T < Tc it fails to tunnel from its early state localized in 
the shallower well to the deeper well as the kinetic en- 
ergy is decreased. This result is qualitatively the same 
for both real-time and quantum Monte Carlo annealing, 
excepting for the fact that the Tc's are different in the 
two cases. 

The realization of QA consists of employing adjustable 
quantum fluctuations into the problem instead of a ther- 
mal one (Kadowaki and Nishimori 1998, Amara et al 
1993, Finnila et al 1994). In order to do that, one needs 
to introduce an artificial quantum kinetic term T{t)TCkin-i 
which does not commute with the classical Hamiltonian 
Tic representing the cost function. The coefficient T is 
the parameter which controls the quantum fluctuations. 
The total Hamiltonian is thus given by 



kin • 



(15) 



The ground state of Htot is a superposition of the eigen- 
states of He- For a classical Ising Hamiltonian of the 
form (2), the corresponding total quantum Hamiltonian 
might have the form (9), where He = ^'l2i>j''-^ ^i^j 
and Hkin = —J^i'^f- Initially F is kept very high so 
that the Hkin dominates and the ground state is trivially 
a uniform superposition of all the classical conflgurations. 
One starts with that uniform superposition as the initial 
state, and slowly decreases F following some annealing 
schedule, eventually to zero. If the process of decreasing 
is slow enough, the adiabatic theorem of quantum me- 
chanics (Sarandy et al 2004) assures that the system will 
always remain at the instantaneous ground state of the 
evolving Hamiltonian Htot- When F is finally brought to 
zero, Htot will coincide with the original classical Hamil- 
tonian He and the system will be found in the ground 
state of it, as desired. The special class of QA algo- 
rithms, where strictly quasi-stationary or adiabatic evo- 
lutions are employed are also known as Quantum Adia- 
batic Evolution algorithms (Farhi et al 2001). 

Two important questions are how to choose an appro- 
priate Hkin and how slow the evolution needs to be in 
order to assure adiabaticity. According to the adiabatic 
theorem of quantum mechanics, for a non-degenerate 
spectrum with a gap between the ground state and first 
excited state, the adiabatic evolution is assured if the 
evolution time r satisfies the following condition- 
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\(f)o{t)) and Ict^iit)) being respectively the instantaneous 
ground state and the first excited state of the total Hamil- 
tonian Htot, and A(<) the instantaneous gap between 
the ground state and the first excited state energies (see 
Sarandy et al 2004). One may wonder that while en- 
tering the ordered phase (F < Fc) from the para phase 
(F > Fc) in course of annealing, the gap A may vanish 
at the phase boundary (F = Fc) in the N ^ oo limit. In 
fact, in such a case, QA cannot help anyway in finding 
the ground state of an infinite system. However, for any 
finite sample, this gap is very unlikely to vanish for a 
random system, and QA may still work out nicely. 

However, it is impossible to follow, even for a finite 
N, the evolution of a full wave function in a classical 
computer using polynomial resources in general, since 
it requires tracking the amplitudes of all the basis vec- 
tors (all possible classical configurations), whose number 
grows exponentially with system-size N. Such an adi- 
abatic evolution may be realized within polynomial re- 
sources only if one can employ a quantum mechanical 
system itself to mimic the dynamics. However, one may 
employ quantum Monte Carlo methods to simulate some 
dynamics (not the real-time quantum dynamics) to sam- 
ple the ground state (or a mixed state at low enough 
temperature) for a given set of parameter values of the 
Hamiltonian. Annealing is done by reducing the strength 
F of the quantum kinetic term in the Hamiltonian from a 
very high value to zero following some annealing schedule 
in course of simulation. In case of such Monte Carlo an- 
nealing algorithm, there is no general bound on success 
time T like the one provided by the adiabatic theorem 
for true Schrodinger evolution annealing. Here we will 
separately discuss the results of real-time QA and Monte 
Carlo QA. Apart from these quasi-stationary quantum 
annealing strategies, where the system always stays close 
to some stationary state (or low-temperature equilibrium 
state), there may be cases where quantum scatterings 
(with tunable amplitudes) are employed to anneal the 
system (Das et al 2005). 



A. Quantum Monte Carlo Annealing 

In quantum Monte Carlo annealing, one may employ ei- 
ther a finite (but low) temperature algorithm, or a zero- 
temperature algorithm. Most of the Monte Carlo QA 
(Das and Chakrabarti 2005, Santoro and Tosatti 2006) 
are done using a finite temperature Monte Carlo, namely, 
the Path Integral Monte Carlo (PIMC), since its imple- 
mentation is somewhat simpler than that of other zero- 
temperature Monte Carlo methods . 
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Among the other zero-temperature Monte Carlo 
methods used for anneahng are the zero-temperature 
transfer-matrix Monte Carlo (The chapter by Das and 
Chakrabarti in Das and Chakrabarti, 2005) and the 
Green's function Monte Carlo (Santoro and Tosatti 
2006). However, these algorithms suffer severely from dif- 
ferent drawbacks, which renders them much slower than 
PIMC algorithms in practice. 

The Green's function Monte Carlo effectively simulates 
the real-time evolution of the wave function during an- 
nealing. But to perform sensibly, it often requires a guid- 
ance that depends on a priori knowledge of the wave 
function. Without this guidance it may fail miserably 
(Santoro and Tosatti 2006). But such a priori knowledge 
is very unlikely to be available in the case of random op- 
timization problems, and hence so far the scope for this 
algorithm seems to be very restricted. 

The zero-temperature transfer-matrix Monte Carlo 
method, on the other hand, samples the ground state 
of the instantaneous Hamiltonian (specified by the given 
value of the parameters at that instant) using a projec- 
tive method, where the Hamiltonian matrix itself (a suit- 
able linear transformation of the Hamiltonian that con- 
verts into a positive matrix, in practice) is viewed as the 
transfer-matrix of a finite-temperature classical system of 
one higher dimension (Das and Chakrabarti 2005). But 
the sparsity of the Hamiltonian matrix for systems with 
local kinetic terms leaves the classical system highly con- 
strained and thus very difficult to simulate efficiently for 
large system sizes. 

The PIMC has so far been mostly used for QA. The 
basic idea of path-integral Monte Carlo rests on Suzuki- 
Trotter formalism (see Appendix 1), which maps the par- 
tition function of a rf-dimensional quantum Hamiltonian 
Ti to that of an effective classical Hamiltonian Tieff in 
{d-\- l)-dimension. Quantum annealing of a Hamiltonian 
TLiot using PIMC consists of mapping Tiioi to its equiva- 
lent classical one and simulate it at some fixed low tem- 
perature so that thermal fluctuations are low. The quan- 
tum fluctuations are reduced from some very high initial 
value to zero finally through reduction of V in course of 
simulation. Clearly, the simulation dynamics is not the 
true Schrodinger evolution of the system, and also it can- 
not simultaneously see the whole configuration space as 
does a dclocalizcd wave hmction. It feels the landscape 
locally and makes moves just like a classical system. The 
first attempt of quantum annealing using PIMC was done 
by Kadowaki and Nishimori (1998) for solving TSP and 
an extensive use of the technique to explore multitude 
of problems has been made by the group of Santoro and 
Tosatti (2006). Here we will discuss some PIMC quan- 
tum annealing results briefly for few different systems. 



1. A Short-Range Spin Glass 

Quantum annealing of an EA spin glass in 2-dimension 
(square lattice) using transverse field (see Eq. (9)) for 



large lattice size (up to 80 x 80) using PIMC turns out 
to be much more efficient compared to thermal anneal- 
ing (CA) of the same system in finding the approxi- 
mate ground state (Santoro et al 2002, Martohak et al 
2002). The quantity which is measured is the resid- 
ual energy €res{T) = E{t) — Eq, Eq being the true 
ground state energy of the finite system (calculated by 
the Spin Glass Server using the so called Branch-and- 
Cut algorithms; see Iittp : // www.informatik.uni — 
koeln.de/lsjUenger/research/spinglass/) and E{t) be- 
ing the final energy of the system after reducing the trans- 
verse field strength F linearly within time r from some 
suitably large initial value to zero. Here r is a fictitious 
time given by the number of Monte Carlo steps. 

Classically, for a large class of frustrated disordered 
system it can be shown, using very general arguments 
that the residual energy decreases following some power 
law in the logarithm of the annealing time r, namely, 
Eres ~ (logr)^'', with C < 2 (Huse and Fisher 1986). In 
PIMC, the partition function of a d-dimensional quantum 
system is mapped to an equivalent d-|-l-dimensional clas- 
sical system (known as Suzuki- Trotter mapping). The 
effective d + 1-dimensional system, is obtained by repli- 
cating the classical part of the original system (say, the 
Ising interaction part of the transverse Ising system with- 
out the transverse field term) with all its interactions 
(including the disorders) intact along the extra higher 
dimension. The coupling between the spin in different 
replica depends on the quantum kinetic term in the orig- 
inal d-dimensional system. 

However, the PIMC annealing results show that the 
quantum effect (taken into account through Suzuki- 
Trotter mapping) does not change the basic relaxation 
behavior of eresij). But still a dramatic improvement 
in evaluation time is achieved since it turns out that the 
value of the exponent ( can be much higher {( = 6) for 
QA than the Huse-Fisher bound of < 2 for classical 
annealing. This is a tremendous improvement in compu- 
tational time if one thinks in terms of the changes in r re- 
quired to improve e,.,.,, equally by some appreciable factor 
in the respective cases of classical and quantum anneal- 
ing. An interesting asymptotic comparison for results of 
QA and CA for an 80 x 80 lattice shows that to reach a 
certain value of eres, the PIMC-QA would take one day 
of CPU time (for the computer used) whereas CA would 
take about 30 yrs (Santoro et al 2002, Martohak et al 
2002). The result would not be much different in this case 
if the real-time Schrodinger evolution would have been 
followed, as has been argued using the Landau-Zener cas- 
cade tunneling picture (Santoro et al 2002, Martonak et 
al 2002). 

Landau-Zener tunneling theory gives an estimate of the 
probability of tunneling non-adiabatically from a lower 

to a higher level when the system encounters an avoided 
level crossing between the levels during time evolution. 
Let the gap between two energy levels \a{t)) and \b{t)) 
of a time-dependent Hamiltonian vary linearly with time 
(gap A = at) and encounter an avoided level crossing. 
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Here the levels are energy levels of the classical part of 
the Hamiltonian (say, the potential energy levels of an 
Ising model). Let the system be made to evolve in such a 
way that the characteristic time it spends while passing 
through the crossing region is t. Let there be a quantum 
tunneling field F that induces transitions between the 
levels. Then if the system evolved is at the lower branch 
I a) before encountering the avoided level crossing, the 
probability that it tunnels to the higher branch |6) while 
passing through the crossing decreases with the time t 
as P{t) = exp(r/rr), where rp = {haV) / {2tt A^^^J , 
Amin being the minimum value of the gap attained at 
the avoided level crossing. For spin glass-like systems 
with non-zero gap, treating multiple level-crossings, each 
with small A, as a cascade of independent Landau-Zener 
tunneling, one may argue that residual energy goes as 
Eres ^ (logT)^^Q, whcrc (q IS essentially greater than 
the bound C 5: 2 for thermal annealing, and might be as 
high as 6 (Santoro et al 2002, Martohak et al 2002). 



2. The Traveling Salesman Problem 

Quantum annealing of TSP with a random metric (i.e., 
the distance dij between the i-th and the j-th city is cho- 
sen randomly) in infinite dimension using PIMC was also 
found to be more efficient than CA in finding an approx- 
imately minimal tour (Martohak et al 2004). An iV-city 
tour in a random TSP problem can be represented by a 
configuration of N'^ constrained Ising spins, and the tour 
length, i.e., the He, by the Hamiltonian (8). Now to do 
the annealing, one needs to introduce a set of moves (spin 
flip operations) that satisfies the constraints. Classically 
one very important class of moves is the 2-opt moves, 
which starting from a valid tour, can generate all possi- 
ble tours without generating any invalid link. Let a valid 
tour contains two links i — s- j and fc — > Z. A 2-opt move 
may consist of removing those two links and establishing 
the following links: i — > fc and j — > I (here i denotes i- 
th city). Classical annealing of the Hamiltonian can be 
done by restricting the Monte Carlo moves within 2-opt 
family only. However, for the quantum case, one needs 
to design a special transverse field (non-commuting spin- 
flip term) to enforce the constraints. It can be realized 
by a spin-flip term of the form S^/. S^^ S^j S^^ ^.^ , 
where the operator S~j flips down (-1-1 —1) the Ising 
spins S^j and S^^ when they are in -f 1 state, and sim- 
ilarly for the flip-up operators S^y However, to avoid 
the Suzuki- Trotter mapping with these complicated ki- 
netic terms, a relatively simple kinetic term of the form 
■Hmu = -r(t) E<^.J> {St^,J> + H.C) is used for the 
quantum to classical mapping. The Monte Carlo moves 
were kept restricted within the 2-opt family to avoid in- 
valid tours. The results were tested on an instance of 
printed circuit board with N = 1002. PIMC-QA was 
seen to do better than the CA and also much better than 
standard Lin-Kernighan algorithm, as shown in Fig. 3. 
The relaxation behavior of residual path-length for TSP 
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FIG. 3 Average residual excess length found after CA and QA 
for a total time r (in MC steps), for the N = 1002 instance prl002 
of the TSPLIB.The dashed horizontal line represents the best out 
of 1000 runs of the Lin-Kernighan algorithm. QA is clearly faster 
than CA (taken from Battaglia et al 2005). 



(see Fig. 3) is found to be quite similar to that of the 
2-d EA glass discussed earlier (see Battaglia et al 2005). 
This indicates that a random TSP also has spin glass-like 
potential energy landscape (PEL), as has already been 
hinted by the replica analysis study of the problem dis- 
cussed in an earlier section (Sec. IIC.2). In that case the 
Landau-Zener tunneling picture would also be applica- 
ble for TSP, and little improvement can be expected by 
following the real-time Schrodinger dynamics instead of 
Monte Carlo methods. 

There are cases where PIMC-QA is not as successful 
as CA. An example (discussed below) of such a prob- 
lem belongs to the class of K- Satisfiability problem (or 
K-SAT). In a /-iT-SAT problem, there is a given Boolean 
function, which is the Boolean sum (connected by logical 
OR operations) of a number of clauses, each of which is a 
Boolean product (connected by logical AND operations) 
of K Boolean variables (binary variables taking values 
or 1 only) taken random from a given set of p variables 
(same variable may occur in various clauses simultane- 
ously). Given such a function, the task is to find an as- 
signment for the Boolean variables for which the number 
of violated clauses [the clauses assigned with the value] 
is the the minimum. The studies on this class of problems 
is extensive and it even includes the connection between 
its hardness and satisfiability-unsatisfiability phase tran- 
sition (Monasson et al 1999). Remarkable progress has 
been made in formulating faster algorithms for solving 
it (see Mezard et al 1987, Hartmann and Rieger 2002) 
based on these understandings. 

But in the case of the random 3-SAT problem (K- 
SAT problem with K= 3) using linear schedule for de- 
creasing F, PIMC-QA gives much worse results than CA 
(Battagha et al 2004). Both CA and PIMC-QA are worse 
than ad- hoc local search heuristics like WALKS AT. In 
the case of application of PIMC-QA in image restora- 
tion, on the other hand, the performance is exactly the 
same as that of CA (Inoue 2005). Even for a particle in 
a simple double-well potential, it seems that with naive 
Monte Carlo moves, PIMC-QA can produce results which 
are much worse than that one could obtain from real-time 
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Schrodinger evolution of the system. There is in fact no 
general prescription to choose the right moves that will 
do the job, unless one has a precise idea about the PEL 
of the problem. The choice of the kinetic term to be 
introduced into the problem is somewhat arbitrary but 
the performance of the algorithm depends on crucially on 
that. It has been found that a relativistic kinetic term 
can do a better job than a non-relativistic one in the case 
of a particle in a double-well potential (Battaglia et al 
2005). PIMC-QA also suffers from difficulty in calculat- 
ing the Suzuki- Trotter equivalent of the quantum Hamil- 
tonian with arbitrary kinetic term designed to satisfy the 
constraints of the problem (Marton'ak et al, 2004). The 
constraints may be taken care of while making Monte 
Carlo moves, but that may not always produce expected 
results. Finally, presence of finite temperature in the 
problem does not allow one to focus exclusively on the 
role of the quantum fluctuations in the problem. Above 
all, like any other Monte Carlo method, PIMC-QA is 
going to do worse if the number of reasonably good ap- 
proximate solutions are only few, and there is no overall 
gradient in the landscape to guide towards them. 

3. Random Field Ising Model: How a Choice of Kinetic Term 
Improves Annealing Results 

QA algorithms enjoy an extra flexibility which a ther- 
mal annealing algorithm cannot. A QA algorithm can 
have a host of choices for its kinetic term. A good choice 
can bring about a lot of betterments. This point is illus- 
trated nicely by Morita and Nishimori (see Morita and 
Nishimori 2007a) for QA of random field Ising model, by 
introducing a ferro-magnetic transverse field interaction, 
in addition to the conventional single-spin-flip transverse 
field term (as present in the Hamiltonian (9)). 

The Hamiltonian of the random-field Ising model with 
the standard single-spin-flip transverse term is given by 

n{t) = nc + n^KL (is) 

where 

JV JV 

nc = -JY.sts'^-Y.htst, (19) 

<ij> i=l 

ftf being the on site random field assuming values -|-1 or 
— 1 with equal probabilities and < ij > denotes sum over 
the nearest neighbor on a 2-dimensional square lattice, 
and 

n 

w^L = ^(^)E^^ (20) 

i=l 

The result of QA in such a system in not satisfactory 
when J is much larger than /if (Sarjala et al 2006). If a 
ferromagnetic transverse term of the form 

JV 

wiii = -m E ^i^l (21) 

<ij> 



is added to the Hamiltonian (18), the result of QA is seen 
to improve considerably (Morita and Nishimori 2007a). 
This happens (as indicated by exact diagonalization re- 
sults on small system) because the ferromagnetic trans- 
verse field term effectively increases the gap between the 
ground state and the first excited state and thus decreases 
the characteristic timescale for the system. This is an ex- 
ample of how one can utilize the flexibility in choosing 
the kinetic term in QA to formulate faster algorithms. 
This also indicates how the knowledge of the phase dia- 
gram of the system, the position of the quantum critical 
point in particular (where the gap tends to vanish), helps 
in choosing additional kinetic terms and thus allows for 
the annealing paths that can avoid, to some extent, the 
regions of very low gap. 



B. Quantum Annealing Using Real-time Adiabatic 
Evolution 

QA is basically the analog version of quantum compu- 
tation. Like the conventional analog quantum computa- 
tion, the hardware realization of adiabatic quantum an- 
nealing is rather problem specific. But once realized, it 
follows the real-time Schrodinger evolution, whose exact 
simulation is always intractable (run time grows expo- 
nentially with the system size) for classical computers 
and often also even for digital quantum computers (see 
Nielsen and Chuang 2000). The annealing behavior with 
the real-time Schrodinger evolution is hence an important 
issue and may show features distinctly different from any 
Monte Carlo annealing discussed so far. 

The first analog algorithm in this line was due to Farhi 
and Gutmann (1998) for solving Grover's search problem. 
Grover (1997) showed that quantum mechanical search 
can reduce an 0{N) classical search time to 0{^/N) time 
in finding out a marked item form an unstructured data- 
base. In the analog version, the problem was to use 
quantum evolution to find a given marked state among 
N orthonormal states. The algorithm was formulated 
in the following way. There were N mutually orthogo- 
nal normalized states, the i-th state being denoted by 
Among all of them, only one, say, the u;-th one, has 
energy E ^ and the rest all have zero energy. Thus 
the state \w) is "marked" energetically and the system 
can distinguish it thus. Now the question is how fast the 
system can evolve under a certain Hamiltonian (our or- 
acle) so that starting from an equal superposition of the 
N states one reaches It was shown (Farhi and Gut- 
mann 1998) that by evolving the system under a time- 
independent Hamiltonian of the form 

ntot = E\w){w\+E\s){s\, (22) 

where 

I^) = ;4EN)' (23) 
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no improvement over Grover's ^/N speedup is possible. 
Later the problem was recast in the form of a spatial 
search (Childs and Goldstone 2004), where there is a d- 
dimensional lattice and the basis state \i) is localized at 
the i-th lattice site. Similarly, as before, the on-site po- 
tential energy E is zero everywhere except at \w). where 
it is 1. The objective is to reach the marked state start- 
ing from the equal superposition of all the |i)s. The ki- 
netic term is formulated through the Laplacian of the lat- 
tice, which effectively introduces uniform hopping to all 
nearest neighbors from any given lattice site and is kept 
constant. The model is in essence an Anderson model 
(Santoro and Tosatti, 2006) with only a single-site disor- 
der of strength 0(1). Grover's speed up can be achieved 
for d > 4 with such a Hamiltonian and no further bet- 
terment is possible. The algorithms succeeds only near 
the critical value of the kinetic term (i.e., at the quantum 
phase transition point). An adiabatic quantum evolution 
algorithm for Grover search was formulated in terms of 
an orthonormal complete set of /-bit Ising-like basis vec- 
tors, where the potential energy of a given basis vector 
(among the 2' ones) is 1, and for rest of all, it is 0. The 
kinetic term is just the sum of all single bit-flip terms 
(as the transverse-field term in Eq. (9)). It connects 
each basis vector to all those that can be reached from it 
by a single bit flip. The kinetic term is reduced from a 
high value to zero following a linear schedule. A detailed 
analysis in light of the adiabatic theorem showed that 
one cannot even retrieve Grover's speedup sticking to 
the global adiabatic condition with fixed evolution rate as 
given by Eq. (16). This is so because the minimum value 
of the gap goes as ^min ~ 2.2~'/^ (Farhi et al 2000), 
and the spin-flip kinetic term being local, the numera- 
tor \(Htot)\ of the adiabatic factor (see Eq. (16)) is at 
best 0{l). However, the Grover's speedup can be recov- 
ered if the condition of adiabaticity is maintained locally 
at every instant of the evolution and the rate is accel- 
erated accordingly whenever possible (Roland and Cerf, 
2001). Adiabatic QA following real-time Schrodinger 
evolution for satisfiability problems also gives strikingly 
different result (for smaller system-size, though) com- 
pared to that obtained using PIMC-QA. Adiabatic QA of 
an NP-complete problem, namely, the exact cover prob- 
lem (as described below) is studied for small systems, 
where a quadratic system-size dependence was observed 
(Farhi et al 2001). In this problem also, the basis vectors 
are the complete set of 2' orthonormal /-bit basis vectors, 
denoted by {|2;i)||2:2)...|2;;)}, each Zi is either or 1. The 
problem consists of a cost function Tic which is the sum 
of many 3-bit clause functions hciizi, Zj, Zk), each act- 
ing on arbitrarily chosen bits Zi, Zj and Zk- The clause 
function is such that hci{zi, Zj,Zk) =0 if the clause CI is 
satisfied by the states \zi), \zj) and \zk) of the three bits, 
or else hci{zi, zj, Zk) = 1. The cost Hamiltonian is given 
by Ti-c = J2ci ^ci- Thus if a basis state \{zi}) dissatis- 
fies p clauses, then 'Hc\{zi\) = p|{zi}). The question is 
whether there exists a basis vector that satisfies all the 
clauses for a given He- There may be many basis vec- 



tors satisfying a clause. All of them will be the ground 
state of He with zero eigenvalue. If the ground state has 
a non-zero (must be a positive integer then) eigenvalue, 
then it represents the basis with lowest number of vio- 
lated clauses, the number being given by the eigenvalue 
itself. The total Hamiltonian is given by 

■Htot{t) = (l - 7) '^k^n + ^-Uc, (24) 

where T-Lkin is again the sum of all the single bit-flip op- 
erators. The initial state at f = is taken to be the 
ground state of Hkini which is an equal superposition of 
all basis vectors. The system is then evolved according 
to the time-dependent Schrodinger equation up to f = r. 
The value of r required to achieve a pre-assigncd success 
probability are noted for different system sizes. The re- 
sult showed a smooth quadratic system size dependence 
for / < 20 (Farhi et al 2001). The result is quite encourag- 
ing (since it seems to give a polynomial time solution for 
an NP-complete problem for small system size), but does 
not really assure a quadratic behavior in the asymptotic 
(/ ^ oo) limit. 

A quadratic relaxation behavior (e^es ~ l/r^, ~ 2) 
was also reported in (Suzuki and Okada 2005a, 2005b) 
for real-time adiabatic QA (employing exact method for 
small systems and DMRG technique for larger systems) 
of a one-dimensional tight-bonding model with random 
site energies and also for a random field Ising model on 
2-d square lattice. In the first case, the kinetic term 
was due to hopping between nearest neighbors, while in 
the second case it was simply the sum of single-spin flip 
operators as given in Eq. (9). Wave function annealing 
results using similar DMRG technique for spin glass on 
ladder has been reported by Rodriguez-Laguna (2007). 

It has been demonstrated that for finite-ranged sys- 
tems, where the interaction energy can be written as a 
sum of interaction energies involving few variables, quan- 
tum adiabatic annealing and thermal annealing may not 
differ much in efficiency. However, for problems with 
spiky (very high but narrow) barriers in the PEL (which 
must include infinite-range terms in the Hamiltonian, 
since in finite-range systems, barrier heights can grow 
at best linearly with barrier width), quantum annealing 
does much better than CA (Farhi and Goldstone 2002) as 
has been argued by Ray et al (Ray et al 1989). Quantum 
adiabatic evolution has been employed earlier to study 
the SK spin glass in a transverse field; stationary charac- 
teristics over a range of the transverse field are calculated 
by varying the transverse field adiabatically (Lancaster 
and Ritort 1997). 

As pointed out earlier in the case of PIMC-QA, the 
kinetic term plays a crucial role in adiabatic QA also. 
This is mainly because the gap A(r) between the ground 
state and first excited state depends on the form of Tikin ■ 
For certain choice of Tikin i it might become very small for 
some r. One then can switch on some additional kinetic 
term (adjusted by another parameter F') to increase the 
dimensionality of the phase diagram and avoid the small 
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gap zone by taking a suitable by-pass. In fact such a 
strategy can convert a failure into success (Farhi et al 
2002). 

So far it has been argued that the advantage of quan- 
tum annealing derives mainly from the fact that quantum 
tunneling can penetrate through very high but narrow 
barriers, which are very hard to jump over thermally. In 
this way the cost-configuration landscapes more acces- 
sible to local moves. This is a key feature that works 
even in the case of quantum Monte Carlo methods like 
PIMC-QA, where one finally samples only a very small 
section of the configuration space. However, another re- 
markable advantage that quantum mechanics provides is 
the ability to "sense" the whole configuration space si- 
multaneously through a delocalized wave function. This 
sensing is largely handicapped by presence of random dis- 
order in the system, because the wave function tends to 
localize in many places, often being unable to to pick up 
the deepest well very distinctly. But this feature may 
be utilized in searching a golf-course- like PEL, namely a 
flat landscape with very deep and narrow wells occurring 
very rarely. If there is only one deep well, the problem is 
just the spatial version of the Grover's problem. As we 
have seen, if the depth x of the well is C(l) (indepen- 
dent of N), then no more than 0{\fN) speed-up can be 
achieved. This can be interpreted as the inability of an 
exceedingly large system to sense a given small wound. 
So one might want to make the depth of the well large 
enough so that it cannot be scaled away as the system- 
size is increased. 

We first consider the case where the well-depth goes as 
X ~ aN , where a is some constant. For a well depth of 
such order for a spatial search at infinite dimension (ki- 
netic term is infinite-ranged) one gets a stunning result: 
The evolution time that guarantees any given probability 
of success becomes independent of N (Das 2007; see Fig. 
4). However, if x scales like x ~ ^ with 7 > 1, then 
the speed-up is lost again - a consequence of quantum- 
mechanical non-adiabaticity. We formulate the problem 
in the following way. As in the case of a spatial search 
discussed above, we denote the state localized at the i- 
th cite by \i) . All the cites except \w) has zero on-site 
potential energy, while |w) have an on-site potential well 
of depth x(0- The system is embedded in infinite di- 
mensions and thus there is a uniform tunneling term F 
between any two sites. To do quantum annealing, we 
evolve x(t) from zero to its final value xo keeping F fixed 
at some moderate constant value. The total Hamiltonian 
is given by 
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FIG. 4 Panel (a) in the figure shows numerical verification of the 
A'-independence of the minimum time Tmin to achieve success prob- 
ability P{\w)) = 0.33. Initially, the system is delocalized equally 
over all sites and evolves with time according to time dependent 
Schrodinger equation with the Hamiltonian (25). As expected from 
exact analytical result for the adiabaticity condition, it is seen that 
Tmin becomes independent of N . Panel (b) shows the variation 
of final probability P{\w)) of finding the state \w) with annealing 
time T for different final value of F, for N = 10^. 



which means that the run time becomes independent of 
N , and the ground state has 0(1) overlap with the tar- 
get state We also confirm it numerically by evalu- 
ating r,„i„, which is the minimal r required for obtain- 
ing a success probability P{\'w)) — 0.33 for different N 
through many decades. Here we have chosen a moder- 
ate constant F, and have evolved the well depth x with 
time. The evolution is computed by solving the time- 
dependent Schrodinger equation numerically and Tmin is 
obtained with an accuracy of 10^^ by employing a bi- 
section scheme. The results (Fig. 4a) clearly show that 
P{\w)) tends to becomes independent of as becomes 
larger and larger. This is completely in accordance with 
the analytical result (see Eq. (26)). 

The relaxation behavior for large A^ for a given anneal- 
ing time r depends on the value of F (see Fig. 4b). If 
F is too small, the system takes a longer time to feel the 
changes in the landscape, and hence the adiabatic relax- 
ation requires a longer time (the adiabatic factor becomes 
bigger; see Eq. (26)). On the other hand, if F is too large, 
the ground state itself is pretty delocalized, and hence the 
final state, though closer to the ground state, has again 
a small overlap with the target state \w). For F = 0.5, 
the result is best. For higher and lower values of F, the 
results are worse, as shown in Fig. 4(b). The relaxation 
behavior is seen to be linear with the annealing time r. 



ntot = ~XQ[^-t/T)\w){w\ 



|i)(j|. (25) C. Annealing of a Kinetically Constrained System 



The time-dependent eigen problem can be solved for 
the above Hamiltonian. It can be shown exactly that if 
Xo aN, then in the N ^ 00 limit, the adiabatic factor 

\{'Htot)\max OL 



^2 . 

min 



(26) 



The adiabatic theorem of quantum annealing assures 
convergence of a quantum algorithm when one starts with 
the initial (trivial) ground state of the Hamiltonian and 
evolves slowly enough so that the system is always in the 
ground state of the instantaneous Hamiltonian. How- 
ever, the benefit of tunneling may be extended even in 
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cases where one does not precisely know the eigenstate of 
the initial Hamiltonian (say for a given unknown PEL) 
and hence is unable to start with it. One might rather 
start with a wave-packet (a superposition of many eigen- 
states of the Hamiltonian) that explores the potential en- 
ergy landscape. Quantum tunneling will still allow it to 
move more easily through the PEL than a classical par- 
ticle if the landscape has many high but narrow barriers. 
A semi-classieal treatment of such a non-stationary an- 
nealing has been discussed in the context of a kinetically 
constrained system (KCS) (Das et al 2005). 

We demonstrate here the effectiveness of quantum an- 
nealing in the context of a certain generalized Kinetically 
Constrained Systems (KCS) (Predrickson and Andersen 
et al 1984). KCSs are simple model systems having trivial 
ground state structures and static properties, but a com- 
plex relaxation behavior due to some explicit constraints 
introduced in the dynamics of the system. These sys- 
tems are very important in understanding how much of 
the slow and complex relaxation behavior of a glass can 
be attributed to its constrained dynamics alone, leaving 
aside any complexity of its energy landscape structure. 
In KCSs one can view the constraints to be represented 
by infinitely high energy barriers appearing dynamically. 

We discuss quantum annealing in the context of a KCS, 
which can be represented by a generalized version of the 
East model (Jackie and Eisinger 1991); a one dimen- 
sional KCS. We also compare the results with that of 
thermal annealing done in the same system. The orig- 
inal East model is basically a one-dimensional chain of 
non-interacting classical Ising ('up-down') spins in a lon- 
gitudinal field h, say, in downward direction. The ground 
state of such a system is, trivially, all spins down. A ki- 
netic constraint is introduced in the model by putting 
the restriction that the i-th spin cannot flip if the (« — 1)- 
th spin is down. Such a kinetic constraint essentially 
changes the topology of the configuration space, since the 
shortest path between any two configurations differing by 
one or more forbidden flips is increased in a complicated 
manner owing to the blockage of the 'straight' path con- 
sisting of direct flips of the dissimilar spins. Further, 
the constraint becomes more limiting as more spins turn 
down, as happens in the late approach to equilibrium. 
As a result, the relaxation processes have to follow more 
complex and lengthier paths, giving rise to exponentially 
large timescale (^ e^/-^ ; Jackie and Eisinger 1991). 

In the Das model (Das et al 2005) there is a chain 
of asymmetric double-wells (each with infinite boundary 
walls), each having a particle localized within them. The 
asymmetry is due to an energy difference of 2/i between 
the two wells of a double well. The particle in one of 
the two (asymmetric) wells can change its location to 
the other well stochastically, either due to the thermal 
fluctuations or the quantum fluctuations present in the 
system. The generalized kinetic constraint is introduced 
by assuming that if the particle in the (i — l)-th double- 
well resides in the lower one of the two wells, then there 
appears a barrier of height x ^-nd width a between the 



two wells of the z-th double-well. In such a situation the 
particle in the i-th double- well has to cross the barrier in 
order to change its location from one well to the other. 
On the other hand, if the particle of the {i — l)-th is 
in its upper well, there is no such barrier to cross for 
the particle to go from one well to the other. Following 
the approximate mapping done in case of a symmetric 
double- well (Chakrabarti et al 1996), this model can be 
approximately represented by a generalized version of the 
East model, where each Ising spin is in a local longitu- 
dinal field h in the downward direction. The spin at the 
i-th site sees a barrier of height x and width a between its 
two energy states when the (i — l)-th spin is down, while 
no such barrier occurs for the z-th spin when the (i — 1)- 
th spin in up. This kinetic constraint is the same in both 
cases irrespective of whether the dynamics is classical or 
quantum. 

When the dynamics of the particle is due to quantum 
fluctuations, the tunneling probabilities come from the 
following semi-classical picture of scattering of a particle 
in a double- well with inflnitely remote outer boundaries. 
If a particle is put in one of the wells of such a double- 
well with some kinetic energy (actually the expectation 
value) r, then it will eventually be scattered by the sep- 
arator (a barrier or step) between the two wells. In such 
a scattering, there is a finite probability P that the par- 
ticle manages to go to the other well. We calculate P 
from the simple picture of scatterings of a particle by 
one dimensional potentials as described below. In the 
thermal case we take simple Boltzmann probabilities for 
crossing the same barriers. The minimum of the energy 
of the Ising chain (equivalent to the potential energy of 
the chain of the double- wells) trivially corresponds to the 
state with all the spins down, i.e., aligned along the longi- 
tudinal field h (where all the particles are in their respec- 
tive lower wells) . To reach the ground state in the quan- 
tum case, we start with a very large initial value of F and 
then reduce it following an exponential schedule given by 
F = Fo exp {—t/Tq). Here t denotes the time and tq sets 
the effective time scale of annealing. At zero-temperature 
the slow spin flip dynamics occurs only due to the tunnel- 
ing (kinetic energy) term F and hence the system ceases 
to have any relaxation dynamics in the limit F ^ 0. The 
barriers are characterized by a height x and a width a, 
the barrier area being (/ = x x a- Similarly, in thermal 
case, we start with a high initial temperature Tq and re- 
duce it eventually following an exponentially decreasing 
temperature schedule given by T = To exp (— f/rc); tc 
being the time constant for the thermal annealing sched- 
ule. Here, when the {i — l)-th spin is down, the flipping 
probability for the i-th spin exp {^x/T))- Otherwise, 
it flips with probability P = 1 if it were in the up state 
and with Boltzmann probability P = exp(— /i/T) if it 
were in the down state. Here in the quantum case, the 
probability of crossing the barrier depends on g so that 
the barrier-width a plays a role, while in the thermal 
case, only x sets the crossing timescale irrespective of a. 

In the simulation (Das et al 2005), N Ising spins 
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FIG. 5 Comparison between classical and quantum annealing for a 
chain of 5 X 10^ spins (for the same initial disordered configuration 
with rrii ~ 10"'^). We show the results for tq = 1.8 X 10^ (for 
quantum) and tq = 10® (for classical) with /i = 1; a lower tq 
would not produce substantial annealing. Starting from the same 
initial values Fq = Tq = 100, (and g = 100 in the quantum case) we 
observe that classical annealing requires about 10^ steps, whereas 
quantum annealing takes about 10* steps for achieving the same 



final order m 



f 



' 0.92. 



{Si = ±1, i — 1,...,N) were taken on a linear chain 
with periodic boundary condition. The initial spin con- 
figuration is taken to be random, so that the magnetiza- 
tion m — i^/N)J2i Si is practically neghgible (m^ « 0). 
We then start with a tunneling field Tq and follow the 
zero-temperature (semi-classical) Monte Carlo scheme as 
mentioned above, using the spin fiip probabilities P's ap- 
propriate for the four cases I-IV. Each complete run over 
the entire lattice is taken as one time unit and as time 
progresses, F is decreased from its initial value Fq ac- 
cording to r = Foe^*/'^Q. The results of thermal and 
quantum annealing are compared in Fig. 5 for the same 
order of initial value and time constant for F and T (the 
barrier height x is taken to be 1000 in both the cases, 
and g was taken to be 100 in the quantum annealing 
case, or equivalently the barrier width a is taken to be of 
the order of 0.1). It is observed that to achieve a similar 
degree of annealing (attaining a certain final magnetiza- 
tion ruf), starting from the same disordered configura- 
tion, one typically requires much smaller tq compared to 
Tc', typically, tc 10'^ x tq for equivalent annealing (for 
similar optimal values of final order mf ^ 0.92). For an- 
nealing with final order nif ~ 1, we find tc ~ 10* x tq. 
This comparison depends on the barrier characteristics 
(the value of g). 



D. Experimental Realization of Quantum Annealing 

Brooke et al (1999) showed experimentally (see also Aep- 
pli and Rosenbaum 2005) that the relaxation behavior in 
reaching deep inside the glass phase in Fig. 6 depends 
on the path chosen. 

A TISG, realized by a sample of LiHoo.44Yo.56F4 in 
laboratory transverse field (F), was taken from a high 
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FIG. 6 Experimental realization of QA and CA in 
LiHoo.44 Y0.56F4 is illustrated on its phase diagram on the tem- 
perature T and transverse field F (measured by the magnitude 
of the external laboratory field in kOe) plane. The material be- 
haves like a conventional ferromagnet in the region labeled FM, 
and shows slow relaxation in the glassy domain wall state labeled 
G. The two paths of relaxations, from an initial point A to a final 
point C on the phase diagram are shown by arrow-headed lines. 
Along the classical path A— >B— >C (dashed arrow) the transverse 
field is not applied until the end, so that the relaxations observed 
are purely thermal. Whereas, along the quantum path A— >D— >C 
(continuous arrow), there is a segment where the temperature is 
small enough and the transverse field is high, so that the fluctua- 
tions are mainly quantum mechanical. Relaxations observed along 
the quantum path are often found to be much faster than those ob- 
served along the classical path at low enough temperature. Taken 
from Aeppli and Rosenbaum (2005). 



temperature paramagnetic to a low temperature glassy 
phase following two separate paths in the F — T plane 
(see Fig. 6). Along the classical path of cooling (CA), 
the transverse field was kept zero throughout, and was 
switched on only after reaching the final temperature. 
The quantum cooling (QA), on the other hand, was done 
in the presence of a high transverse field, which was low- 
ered only on reaching the final temperature. As the sam- 
ple is cooled, spectroscopy of the sample at different tem- 
peratures (both during CA and QA) was done to reveal 
the nature of the distribution of spin relaxation time 
scales. The QA produced states whose relaxation was 
up to 30 times faster than those produced by the CA at 
low temperature (see Aeppli and Rosenbaum 2005). This 
clearly indicates that quantum tunneling is much more 
effective in exploring the configuration space in the glassy 
phase than thermal jumps (as indicated in Fig. 1). An 
experimental realization of quantum adiabatic anneal- 
ing for 3-bit instances of MAXCUT problem using NMR 
technique has also been reported (see Steffen et al 2003). 
Here, the the smoothly varying time dependent Hamil- 
tonian was realized by the technique of quantum simula- 
tion (see Nielsen and Chuang 2000), where a smooth time 
evolution was achieved approximately (Trotter approxi- 
mation) through the application of a series of discrete 
unitary operations. The results indicate existence of an 
optimal run time of the algorithm. 
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IV. CONVERGENCE OF QUANTUM ANNEALING 
ALGORITHMS 

Here we briefly summarize some important recent results 
derived by Morita and Nishimori (2006, 2007b) on the 
convergence of QA algorithms for TIM systems. The 
results are valid for both the Quantum Monte Carlo and 
the real-time Schrodinger evolution versions of QA. 

The Hamiltonian here is the same as given in (9) with 
a time dependence in the transverse field F = T{t). No 
assumption regarding either the nature of the distribu- 
tion of Jij or the spatial dimensionality is required. In 
order to perform QA at temperature T using PIMC, one 
constructs the Suzuki- Trotter equivalent (see Appendix 
1) d+l dimensional classical system of the d dimensional 
quantum system, and the resulting (classical) system is 
simulated using a suitable inhomogeneous Markov Chain. 
The transverse field T(t) is tuned from a high value to 
zero in course of simulation. It can be shown that at the 
end of the simulation the final distribution converges to 
the ground state of the classical part of the Hamiltonian 
irrespective of the initial distribution (strong ergodicity) 
if 



r{t) > AfTtanh" 



1 



it + 2) 



2/RL 



(27) 



where M is the number of Trotter replicas in the d+l 
dimensional Suzuki- Trotter equivalent system. Here R 
and L are constants depending on the system size N , 
the spin-fiip dynamics appointed for the simulation, the 
temperature T, etc. For large t the bound reduces to 



V{t) > MT 



1 



{t + 2)2/-Ri ■ 



(28) 



V. QUANTUM QUENCHING 

Instead of annealing, one can quench a systems by re- 
ducing the transverse field from Ti{> Fc) to F/(< Fc) 
and follow the relaxation dynamics after that. Here Fc 
denotes the quantum critical point for the system. This 
can help to prepare a state of the system with substantial 
order in a very short time. Recently, quantum quench- 
ing dynamics in different systems, in particular for sud- 
den quenching across a quantum critical point, have been 
studied extensively (Scngupta et al 2005, Calabrese and 
Cardy 2006, Das et al 2006). The (exact) results for the 
after-quench dynamics, following a sudden quench of an 
infinite-range Ising model in transverse field (Das et al 
2006, see Appendix 2) is shown in Fig. 7. The results 
indicate that for a purely quantum system (small S) , one 
might reach a non-stationary oscillatory state with a sub- 
stantial value for the long-time average (over the largest 
timescale) of the order the quenching is 

done from F^ > Fc to F/ = Fc/2. Thus one gets a state 
(though non-stationary) which has considerable order, in 
just one step quenching. However, if the quenching is 
carried out with a high but finite rate, one ends up reach- 
ing an excited state with topological defects (see Kibble 
1976, Zurek 1985, Sengupta et al 2008). 




S=50 
S=100 
S=200 
S=500 



It is remarkable that in contrast to the inverse logarith- 
mic decay of temperature required for convergence of 
CA (see Eq. 1), QA requires only a power-law decay 
of the transverse field. In this sense, quantum anneal- 
ing is much faster than the classical annealing for TIM. 
Similar result can be derived for more general form of 
Hamiltonian (see Morita and Nishimori 2006). However, 
the advantage gained here does not change the complex- 
ity class of an NP-complete problem, since RL is of the 
order of N and hence the convergence time is thus still 
exponential in N. 

For real-time Schrodinger dynamics at T = 0, the 
bound for the decay of the transverse field is again of 
the form 



(29) 



where ^ is exponentially small in N for large N. Though 
the dynamics of PIMC-QA and real-time Schrodinger 
evolution are completely different, the power-law bound 
on the annealing schedule is strikingly similar. 



FIG. 7 Plot of the long-time average (see Appendix 2) O = 
i(^tot)'^) ^ ^ function of Tf/J for different S. In the plot the 
solid (blue), dotted (black), dash-dotted (green) and the dashed 
(red) lines represents respectively the results for S = 50, S = 100, 
S = 200 and S = 500 (color online). O peaks around Tf/J = 0.25 
and the peak value decreases with increasing S. For all plots we 
have chosen Fi/J = 2. (The figure is taken from Das et al 2006). 



VI. SUMMARY AND DISCUSSIONS 

Unlike the gate-based quantum computers (see, e.g., Ek- 
ert and Jozsa 1996, Nielsen and Chuang 2000, Galindo 
and Martin-Delgado 2002), the annealing of a physical 
system towards the optimal state (encoded in the ground 
state of the final Hamiltonian) in the classical limit nat- 
urally achieves analog quantum computation. As dis- 
cussed here, the utilization of the quantum mechanical 
tunneling through the classically localized states in an- 
nealing of glasses has opened up this new paradigm for 
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analog quantum computation of hard optimization prob- 
lems through adiabatic reduction of the quantum fluctu- 
ations. 

We reviewed here the recent success in annealing, or 

optimizing the cost functions of complex systems, utiliz- 
ing quantum fluctuations rather than the thermal fluc- 
tuations (see Santoro and Tosatti 2006 for a more tech- 
nical review). As mentioned already, following the early 
indication in Ray ct al (1989) and the pioneering demon- 
strations, theoretically by Amara et al (1993), Finnila ct 
al (1994), Kadowaki and Nishimori (1998), Farhi et al 
(2001) and Santoro et al (2002), and experimentally by 
Brooke et al (1999), the quantum annealing technique 
has now emerged as a successful technique for optimiza- 
tion of complex cost functions. The literature, exploring 
its success and also its limitations, is also considerably 
developed at present. 

These are introduced here through discussions of the 
mapping of such hard problems to classical spin glass 
problems, discussions on quantum spin glasses, and con- 
sequent annealing. The physics of classical spin glasses 
(see Sec. II) offers us the knowledge of the of the energy 
or the thermodynamic potentials and their landscape 
structures. Mapping of computationally hard problems 
like the TSP etc problems to their corresponding classical 
spin glass models also helped understanding their com- 
plexity (Sec. IIC). The timescale for tunneling through 
an energy barrier quantum mechanically, involves not 
only the height of the barrier, but also its width; the 
narrower the barrier, the faster the tunneling (for con- 
stant height). Thermal fluctuations, on the other hand, 
see only the barrier height for crossing it. This reduction 
of tunneling time with the barrier width helps relaxing a 
quantum system much faster in some cases. This leads to 
the quantum annealin, a framework for constructing gen- 
eral lieuristics to find approximate solutions of hard opti- 
mization problems. While simulated annealing employs 
the strategy of slow cooling, physical or in simulations, 
to find the ground state of glassy systems, quantum an- 
nealing employs quantum fluctuations (see Sec. III). As 
mentioned before, this effectively reduces Planck's con- 
stant to zero to reach the classically optimized (minimum 
cost) state. This reduction, when done completely adia- 
batically, guarantees that the system will be found in the 
ground state of the classical glass at the end (provided 
there is no crossing of energy levels with the ground state 
in the course of evolution, and one has started initially 
with the ground state of the Hamiltonian) . This has al- 
ready been realized experimentally (see Sec. HID), where 
faster relaxation towards the ground state is achieved by 
reducing the external field (inducing changes in the tun- 
neling field), rather than by reducing the temperature. In 
this way analog quantum computation is realized through 
a novel route. Recently an equivalence between the adia- 
batic QA and standard gate-based quantum computation 
has also been establised (Aharonov et al 2007, Mizel et 
al 2007). 
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VII. APPENDIX 

1. Suzuki-Trotter Formalism 

Here we illustrate this formalism by applying it to a TIM 
having a Hamiltonian 

N 

W = -r ^ 5f - ^ J,, S^S^ = Hkin + He . (30) 

The canonical partition function of H reads 

Z = Tre-(«^*''+«^)/^. (31) 
Now we apply the Trotter formula 

exp{Ai+A2)= lim [exp Ai/Mexp^z/M]^ , (32) 

M^oo 

even when [Ai, A2] ^ 0. On application of this, Z reads 
^ = J™ Y.('^\[e^p{-n,,JMT) X 



exp(-Wc/MT)]^|si) 



(33) 



Here Sj represent the i-th spin configuration of the whole 
system, and the above summation runs over all such pos- 
sible configurations denoted by i. Now we introduce M 
number of identity operators in between the product of 
M exponentials in Z, and have 



M , 

Z = lim Tr TT(S'i,fe...S'jv,/c| exp ( - 

M— >(X) -'•-^ \ 
fc=l ^ 



MT 



exp 



-He 
MT 



\Sl^k+l---SN,k+l), 



and periodic boundary condition would imply Sn+i,p = 
Si^p. Now, 

]^(S'i,fc...S'jv,/c| exp I -j^^Jij^i^j \Si,k+i—SN,k+i) 
fe=i \ hj / 



= exp 



N M 

2^ MT 

i,j=l ft=l 



Si,kSj,k 
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where Si^k = ±1 arc the eigenvalues of operator (see 
Hatano and Suzuki, 2005), and 



M 



]J(S'i,fe...S']v,fc| exp 



fe=i 





■ 2r ■ 




Qsinh 


_MT 


) 







exp 



1 / r \ ^ ^ 

^ ^ i=i fc=i 



|5'i,fe+i"-'S'jv,fe+i) 



giving the effective classical Hamiltonian (14), equivalent 
to the quantum one in (30). In the above equation M 
should be at the order of 1/T {% = 1) for a meaning- 
ful comparison of the interaction in the Trotter direc- 
tion with that in the original Hamiltonian. For T — > 0, 
M oo, and the Hamiltonian represents a system of 
spins in a ((i-l-l)-dimensional lattice, because of the ap- 
pearance of one extra label k for each spin variable. Thus 
corresponding to each single quantum spin variable Si in 
the original Hamiltonian we have an array of M number 
of classical replica spins Sik ■ This new (time-like) dimen- 
sion along which these classical spins are spaced is known 
as Trotter dimension. 



2. Quantum Quenching of a Long Range TIM 



Let us consider a system of spin ^ 

the Hamiltonian H = E,>j SfS^ - T 
easily be rewritten as 



^ objects governed by 
Sf. It can 



n 



(34) 



where Sf, 



QX 



Ei ^i^d ^ constant 

(J/27V)E^('5f)2 = J/8 from H in (34). The above 
Hamiltonian can again be cast into the simplified form 
Ti = h.Stoti where h = Jm,z — Fx, giving the mean field 
(exact in this long-range limit) equation m = (St^t) = 



If •»-('?) 



Jm 



2vTO^' at T = and z, x de- 
note unit vectors along z and x directions respectively. 
This gives to = for F > Fg and to ^ for F < 
Fc = J/2. Since the model is infinite-ranged one, the 
mean field approximation becomes exact and one can 
readily express S in terms of its polar components as 
S = S{sva.6 cos (f), sin 6 sin (f), cos 9), S being the total an- 
gular momentum. One can immediately utilize the clas- 
sical equation of motion ^ = S x h. Considering the 
above equation for the z and x components, we get 



dO ^ . . , del) 
— = 1 sm 6* and — 
dt dt 



^ cos ^ -|- F cot ^ cos ^. (35) 



Here we have S = N/2. If the system is now quenched 

from above its quantum critical point F > Fc, finally 
to a F/ < Fc, then one can write (see Das et al 2006), 



equating the energies of the states with and without any 
order respectively, F/ = -j cos^ 6 + T f sin9 cos 0. Using 
this, one gets from Eq. (35) 



^F^sin^^-[F^-J/4cos2^"]^ 



smf 



f{e). (36) 



has zeros (turning points) at 9i = 
l-4Fj/J|) and 02 = 7r/2. One can there- 
(cos^^) = Af/V, where 

jgyecos^e/f{e) = 4^8r/(j-2r/)/j and 

fg^ d9/ f{6), giving a behavior shown in Fig. 7. 



This 
sin~^ 

fore obtain {{S^^J') 
V 
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